MODULE My_Tools

INTERFACE
  Subroutine tauchen(STDE,RHO,SZ,PZW,PZW0,PZEE,pi,ZZ,M,mu,MEANW0,STDW0)
integer, parameter :: dd = Selected_real_kind(15,307)
integer :: ZZ
Real (Kind=dd), Intent(In) :: STDE,RHO,M,mu
Real (Kind=dd), Intent(Out) :: SZ(ZZ),PZW(ZZ),PZEE(ZZ,ZZ),pi(ZZ),PZW0(ZZ)
Real (Kind=dd) ::  MEANW0,STDW0
END Subroutine tauchen
END INTERFACE

INTERFACE
SUBROUTINE INTEGRALTAUCH(XA,XB,VINT,MU,STD)
integer, parameter :: dd = Selected_real_kind(15,307)
Real (Kind=dd), Intent(In) :: XA,XB,MU,STD
Real (Kind=dd), Intent(Out) :: VINT
END Subroutine INTEGRALTAUCH
END INTERFACE

INTERFACE
Subroutine NORMAL_value(X,MU,STD, NORMAL)
integer, parameter :: dd = Selected_real_kind(15,307)
Real (Kind=dd), Intent(In) :: X,MU,STD
Real (Kind=dd), Intent(Out) :: NORMAL
END Subroutine NORMAL_value
END INTERFACE

INTERFACE
Subroutine discretize_normal(hmin,hmax,nhc,hmean0,hstdv0,h,fh)
integer, parameter :: dd = Selected_real_kind(15,307)	
Integer, INTENT(IN) :: nhc
Real (Kind=dd), INTENT(IN) :: hmin,hmax,hmean0,hstdv0
Real (Kind=dd), INTENT(OUT) :: h(nhc),fh(nhc)
END Subroutine discretize_normal
END INTERFACE

END MODULE My_Tools
!======================================================================================


Subroutine tauchen(STDE,RHO,SZ,PZW,PZW0,PZEE,pi,ZZ,M,mu,MEANW0,STDW0)
!=============================================================
!DISCRETIZE A NORMAL DISTRIBUTION BY USING THE TAUCHEN(1986)
!==============================================================

!	Input: MEANW(mean of PZW(z)~Ln(MEANW,STDW^2));
!            STDW(standard deviation)
!	       Log(z_t+1)=RHO*Log(z_t)+e_t+1; e~N(0,STDE^2).
!	Ouput: SZ(set of technology shocks),PZW(pdf),
!	       PZE(Transition matrix)
!		   pi is a stationary distribution
USE My_Tools, ONLY: INTEGRALTAUCH,NORMAL_value
Implicit NONE

INTEGER :: I,IZ
integer :: ZZ
integer, parameter :: dd = Selected_real_kind(15,307)
Real (Kind=dd), Intent(In) :: STDE,RHO,M,mu
Real (Kind=dd), Intent(Out) :: SZ(ZZ),PZW(ZZ),PZEE(ZZ,ZZ),pi(ZZ),PZW0(ZZ)
Real (Kind=dd) ::  MEANW,MEANW0,STDW,STDW0,step, SUM1 ,x


MEANW=mu


STDW=STDE/sqrt(1-RHO*RHO)
SZ(ZZ)=M*STDW

SZ(1)=-M*STDW
step=2*STDW*M/real(ZZ-1)
DO IZ=2,ZZ-1
   SZ(IZ)=SZ(IZ-1)+step
End do


CALL INTEGRALTAUCH(SZ(1)-30.0,SZ(1)+step/(2.0),PZW0(1), MEANW0, STDW0)
      
CALL INTEGRALTAUCH(SZ(ZZ)-step/(2.0),SZ(ZZ)+30.0,PZW0(ZZ), MEANW0,STDW0)
	
DO IZ=2,ZZ-1
   CALL INTEGRALTAUCH(SZ(IZ)-step/(2.0), SZ(IZ)+step/(2.0),PZW0(IZ),MEANW0,STDW0)
Enddo
PZW0= PZW0/sum(PZW0)

!==================================
!DISCRETIZE LN(z)~N(MEANW,STDW**2)
!==================================
	
CALL INTEGRALTAUCH(SZ(1)-30.0,SZ(1)+step/(2.0),PZW(1), MEANW,STDW)
      
CALL INTEGRALTAUCH(SZ(ZZ)-step/(2.0),SZ(ZZ)+30.0,PZW(ZZ), MEANW,STDW)
	
DO IZ=2,ZZ-1
   CALL INTEGRALTAUCH(SZ(IZ)-step/(2.0), SZ(IZ)+step/(2.0),PZW(IZ),MEANW,STDW)
Enddo
	


x=0.0

DO IZ=1,ZZ
   DO I=2,ZZ-1
      CALL INTEGRALTAUCH( &
             SZ(I)-MEANW*(1.0-RHO)-RHO*SZ(IZ)-step/(2.0), &
             SZ(I)-MEANW*(1.0-RHO)-RHO*SZ(IZ)+step/(2.0), &
             PZEE(IZ,I),x,STDE)
   Enddo
Enddo
	
DO IZ=1,ZZ
   CALL INTEGRALTAUCH(SZ(1)-30.0, SZ(1)-MEANW*(1.0-RHO)-RHO*SZ(IZ)+step/(2.0), &
                    PZEE(IZ,1),x,STDE)
         
   CALL INTEGRALTAUCH(SZ(ZZ)-MEANW*(1.0-RHO)-RHO*SZ(IZ)-step/(2.0), &
                     SZ(ZZ)+30.0,PZEE(IZ,ZZ),x,STDE)
Enddo

!=============
!NORMALIZATION
!=============
SUM1=0.0
DO I=1,ZZ
   SUM1=SUM1+PZW(I)
Enddo
	
DO I=1,ZZ
   PZW(I)=PZW(I)/SUM1
Enddo
	
DO IZ=1,ZZ
   SUM1=0.0
	DO I=1,ZZ
		  SUM1=SUM1+PZEE(IZ,I)
	Enddo
	
	DO I=1,ZZ
		  PZEE(IZ,I)=PZEE(IZ,I)/SUM1
	Enddo
Enddo
	
DO IZ=1,ZZ
   SZ(IZ)=EXP(SZ(IZ))
Enddo
	

pi=1.0/ZZ
Do i=1,150
	pi=Matmul(pi,PZEE)
Enddo



END Subroutine tauchen
!****************************************************************************
!****************************************************************************
	
SUBROUTINE INTEGRALTAUCH(XA,XB,VINT,MU,STD)
!	==============================================================
!	COMPUTATION OF A NUMERICAL INTEGRAL BY A SERIES OF RECTANGLES
!	WITH EQUAL WIDTHS
!	=============================================================

!	NOTE: Calculate area under a given curve by approximating 
!	      the given region with a series of rectangles with equal
!           widths. The heights of the rectangles are the individual
!           ordinates associated with each X value of the curve. 
!           The areas of the rectangles are summed. 
!	Input: N(Nb of points between xa and xb),
!	       XA(lower bound of integral),XB(upper bound of integral)
!	Ouput: VINT(value of integral)
USE My_Tools, ONLY: NORMAL_value
integer, parameter :: dd = Selected_real_kind(15,307)	
INTEGER, PARAMETER :: N=1000
INTEGER I
Real (Kind=dd), Intent(In) :: XA,XB,MU,STD
Real (Kind=dd), Intent(Out) :: VINT
Real (Kind=dd) ::  DX,X,NORMAL
	
DX=(XB-XA)/real(N)
X=XA
Call NORMAL_value(X,MU,STD,NORMAL)
VINT=NORMAL
	
DO I=1,N
   X=X+DX
	Call NORMAL_value(X,MU,STD,NORMAL)
   VINT=VINT+NORMAL
Enddo

VINT=DX*VINT

RETURN
END Subroutine INTEGRALTAUCH

!*****************************************	
Subroutine NORMAL_value(X,MU,STD, NORMAL)
!================================
!NORMAL DISTRIBUTION
!================================
!	Input: X(random variable),MU(mean),STD(standard deviation)
!	Output:NORMAL(probabilities of X)

integer, parameter :: dd = Selected_real_kind(15,307)
Real (Kind=dd), Intent(In) :: X,MU,STD
Real (Kind=dd), Intent(Out) :: NORMAL
Real (Kind=dd) :: TP,PI

PI=3.141592653589790
TP=STD*(2.0*PI)**(0.50)
TP=1.0/TP
NORMAL=(X-MU)/STD
NORMAL=NORMAL**2.0
NORMAL=-0.50*NORMAL
!print *, NORMAL
!pause
NORMAL=EXP(NORMAL)
NORMAL=TP*NORMAL
RETURN
END Subroutine NORMAL_value


!===================================================================================
Subroutine discretize_normal(hmin,hmax,nhc,hmean0,hstdv0,h,fh)
USE My_Tools, ONLY: INTEGRALTAUCH
Implicit NONE
integer, parameter :: dd = Selected_real_kind(15,307)	
Integer, INTENT(IN) :: nhc
Real (Kind=dd), INTENT(IN) :: hmin,hmax,hmean0,hstdv0
Real (Kind=dd), INTENT(OUT) :: h(nhc),fh(nhc)
Integer :: jhc,index
Real (Kind=dd) :: g,fhfem(nhc),hmean0fem,hstdv0fem
!     INPUTS: minimum and upper bound for grid points, mean of log human capital
!             and standard devition of log human capital. 
!     OUTPUT: grid for human capital h
!             initial distribution (age0) for human capital
!     initialize
h(1)=hmin
h(nhc)=hmax
	
!     This program constructs a grid such that 
g=(h(nhc)-h(1))/real(nhc-1)
do 1 jhc=2,nhc-1
	 h(jhc)=h(jhc-1)+g
1  continue
  
CALL INTEGRALTAUCH(h(1)-30.0,h(1)+ g/(2.0),fh(1), hmean0, hstdv0)

CALL INTEGRALTAUCH(h(nhc)-g/(2.0),h(nhc)+30.0, fh(nhc),hmean0,hstdv0)
	
DO jhc=2,nhc-1
   CALL INTEGRALTAUCH(h(jhc)-g/(2.0),h(jhc)+g/(2.0),fh(jhc),hmean0,hstdv0)
Enddo
fh=fh/sum(fh)


end subroutine discretize_normal

